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Introduction 

A considerable amount of interest has recently developed 
in the correction of spatial image distortions for registration 
of multitemporal remote sensor imagery [1,2,33. The problem 
arises in the case where a scene is imaged at two or more times 
under varying sensor states. It is desired that the multiple 
images be geometrically registered so that when converted to 
digital form they can be analyzed as a multidimensional image 
vector using computer techniques. 

One technique for registering two images is to find corres- 
ponding points in the two images and use these point pairs to 
distort one image to match the other. This problem has two 
characteristicjs which make specification of the correction function 
difficult. The first is the fact that the distortion of one 
image with respect to another is highly sensor dependent. Certain 
sensors introduce a great deal of random spatial distortion while 
others are highly stable. An example of a distortion producing 
sensor is the multispectral line scanner carried by a low altitude 
aircraft. Pitch, yaw, and translational movements of the aircraft 
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introduce corresponding distortions in the imagery. A highly 
stable sensor example is a multispectral scanner or camera carried 
in satellite orbit about the earth. Only slight distortions 
are introduced in instruments of this type. 

The second key problem is that identification of matching 
points or checkpoints, as they will be called, is highly data or 
scene dependent. Matching points cannot be found at regular 
intervals either visually or by correlation by using scene context. 
Road intersections or correlated scene features tend to be found 
at random over an area. In images gathered at widely separated 
epochs the scene may have undergone such drastic change that very 
few matching points can be found. This problem does not exist if 
tick marks or reseau grids, as they are properly called, are imaged 
in coincidence with the scene. This is possible for image forming 
sensors such as cameras or vidicons but not for scanners. 

The problem addressed in this paper is that of determining 
the optimum two dimensional approximation function for image 
distortions when the data or checkpoints are unequally spaced. 
Although the distortion function is a two dimensional function 
of two dimensions, it can be separated into two independent one 
dimensional functions of two dimensions. 

The general statement of the problem is; 

Find : 

A x (x,y) = F x (c,p,x,y) 

A y (x,y) = F y (c,p,x,y) 


( 1 ) 
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such that the distortion functions A , A are as "close" as 

x y 

possible to the true distortions of the subject image. Defini- 
tion of the functional form and a meaning of "close" are two 
problems of equal importance. In the above: 

A (x f y) are the estimated distortion values in two 

x t y 

dimensions 

are the approximating functions 
c checkpoints 

p parameters defining F 

A formal statement of the problem requires assumption of some 
functional form for the true distortion over the two dimensional 
space f(x,y) considered and some form for the error p. Then 
the problem can be stated : 

Let f(x,y) be a real valued continuous function of two 
variables on a set R, and let F(A,x,y) be a real valued 
approximating function depending continuously on x,y^ R 
and on parameter A. Given the error function p, determine 
the parameters A* Q such that 

p [F (A* ,x,y) t f(x,y)] < p[F(A,x), f (x) ] { 2 ) 

for all A Q, where Q is the set of all possible parameter sets. 
The choice of an approximating function is difficult since 
no explicit method exists for making such a choice. Only the 
general statement that the more comolex the variations in the func- 
tion are the more complex the approximating function must be can 
be made with certainty. The choice of error function is also 
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equally undefined. The choice of error measure is often based 

on generally favorable characteristics of certain well known 

functions. The error or distance function is commonly called 

a norm and a common class of norms which will be considered 

here is called the L norm. The L norm of the function f(x) 

P P 

is denoted by L (f) and is defined by: 
p b 


L p (f) = 


| i f (x) | dx 
a 


1/P 


p > 1 


( 3 ) 


Then the best approximation to f (x) is obtained when the 
L distance function is minimized: 

p b p 

m ^ n | | F (A,x,y) - f<x,y)| dx (4) 

a 

The solution A* for p = p^ will in general be different 
for p = p^ and the nature of the approximation varies sharply 
as p is varied. Some well known cases are for p = 1 which is 
the minimum sum of the absolute value norm and p - 2 the least 
squares norm. A third widely used norm is called the Tchebycheff 
norm which is simply the maximum error. An optimum approximation 
in the Tchebycheff norm minimizes the maximum error. The or 
least squares norm is generally preferred above all others 
because of its desirable heavier weighting of large errors more 
than small errors and because of its differentiability and its 
relationship to series of orthogonal functions. The norm will 
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be used in the work described here. 

The form of the approximating function remains as the 
key problem in the study. The functions investigated in 
this study are from two classes: 1. Polynomials, and 2. 

Spline Functions. A large body of information exists on 
polynomial approximations and the work reported here is 
based largely on Ralston [ 4 ] and Rice [5]. Spline functions 
or piecewise continuous polynomials have received relatively 
little attention until the early 1960 's and the work here is 
based primarily on Rice [6] and deBoor [7]. Polynomials are 
studied first and are used to define image distortion over the 
entire two-dimensional image space. As distortions became 
more severe the order of the approximating polynomial becomes 
high and solution problems become severe. Spline functions 
are investigated secondly to determine if lower order polynomials 
fitted together can approximate a higher ordcs distortion with 
less computational difficulty. Theory for the one dimensional 
case is first developed and then the two-dimensional theory is 
developed. Algorithms for generating approximating functions are 
described and application of the techniques to description 
of image distortion in aircraft multispectral scanner imagery 
is described. Comparison of results and recommendations are 
presented in conclusion. 
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One-Dimensional Polynomial Approximation 

Approximation of a function defined by a discrete set 
of points {f i**i , , . . ,n) defined at points using a set of 
polynomials of largest order m {$ j (x) I j“0 / • • . is expressed 

as: m 

F{A,x)* 2 a . 6 . (x) (5) 

j=0 33 

The least squares approximation to {f.j_} is defined by 

the set of coefficients A* such that 
n , 

e A “iJ ^ (fi"F (A,x) ) 2 is a minimum. (u) 

The error can be minimized in the L norm using differentia- 

2 

tion since the I» 2 norm has the desireable property of being 
analytically differentiable. Taking the derivative of 
with respect to each coefficient and setting the result to 
zero gives: 

n m 

•!»_ - - 2 -jSo W=° < 7 > 

9a k 

A unique solution A* is guaranteed since the L norm is a 
strictly convex function of A. This expression creates 
what is known as the normal equations for the least squares 
approximation. Another advantage of the L 2 norm is that the 
resulting normal equations are linear in the parameter space 
A whereas higher order norms result in quadratic, cubic, etc. 
equations in A. The equations can be expressed as: 
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m 


j£ 0 g j k a j" P k k“0, . . ,m 


where : 


and 


n 


g jk“i®i^j (x i^k* x i* j /k»0, , . . ,m 


n 


In matrix form notation this formulation is expressed as: 
G-fg^M^ (x A ) } 


(8) 

(8 a) 
(8b) 

(9) 


(q! . denotes the elements of G not g.. in equation 8a) 
i 3 3 ^ it» 

Then the summation for g^ terms of the matrix G is G G [ 8 J • 

T 

The summation for is expressed as G f, where f is a 
column vector {f^}. Then the expression for the a^ becomes: 


T T 
G A GA«G A f 


( 10 ) 


( 11 ) 


Where A is a column vector of the desired coefficients 
{a^}. The solution for A is obtained by solving the above 
matrix equation for A 

m - Irp 

A*(G G) G f 

The elements of G are the basis functions^ j (x^) evaluated at 
the ordinates {x^} and are independent of f. The elements 
of G are, however, highly dependent on the form of ♦. If 
4> j (x^) is chosen to be x^ then the exponent of the abscissa 
grows as the square of the order of the approximation. The 
matrix G G becomes highly ill-conditioned for values of m 
larger than 5 or 6 and this makes inversion difficult due to 
round-off error in the computational process. To illustrate 

this, assume all values of Xj lie in the interval 0->l. The 

T i 

terms in G G are for $j»x : 


( 12 ) 
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* 


n 


a » E J +k 
9 jk i*l x i 


which is approximately n times a Riemann sum 

* 


g ;Jk “ n /Jx ;j+k dx« j ,k=0, . . . ,m 


( 13 ) 


Which in matrix form is: nH where: 

i 

• • 

m+ 1 


H- 


l i 

ITT 


1 i i 

T T V 


[m+i 2m+i j 

The matrix H is the principal minor of the infinite Hilbert 
matrix and is a classical example of an ill-conditioned 
matrix. An ill-conditoned matrix is one which when its 
largest value is 1 has an inverse with very large values. 

Thus use of the basis functions = should be restricted 
to approximations of low order. 

T 

Note that the elements of G G are cross or inner products 
of the space of basis functions. If the are chosen such 
that 


n 


ii, ^ “ 0 for j/k and 
1 * 3 1 * 1 ft 0 for j-k 


( 14 ) 


then the off diagonal elements of G T G are zero and the 
matrix is easily invertible. Polynomials with this property 
are called orthogonal and are of extreme importance in the 
theory of approximation. Furthermore if 
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for j=k (15 ' 

The <f>^ are called orthonormal and the matrix G T G becomes 
the identity matrix and the inversion problem vanishes. The 
approximation problem becomes 


T 

IA-G 1 f 


<1G) 


which is a simple matrix multiplication. The elements of 
A are n 

A problem arises in that for finite point sets the ortho- 
normal polynomials depend on the number and spacing of the 
points x^. For cases where the points are equally spaced 
the polynomials can be defined parametrically in the number 
of points n and the order of the polynomial m. These poly- 
nomials are known as the Gram polynomials and can be found 
tabulated and their derivation is unnecessary. If, however, 
the points are unequally spaced as is the case for the 
checkpoints defining the image distortion, the Gram poly- 
nomials cannot be used. In this case, a recursive formula 
must be usr*<3 ho generate the polynomials. It can be shown 
by induction (see Ref 4, Page 241) that the following 
relationship generates the orthogonal polynomials: 


♦ j+ l <*) - < x “ c j+i> ♦ j ( x > “ d j * j-i < x) 


(17) 


“10 


Where: 

P Q (x)=l 

P^M-O 

and C , d. are to be determined. The coefficients are 
3 3 

® nd n n 

C k + 1“ i£, 

These polynomials are orthogonal but not normal, thus a 
normalizing factor is. required to compute the coefficients 
for the approximating function. 


m 

F(A,x)= a^j (x) 

( 1 8 ) 

V i J, f iW 



And the normalizing factor is 

v X 

These are in fact the diagonal elements of the matrix G G. 

If the order of the approximating polynomial is low {4 or 5) 
the non-orthogonal normal equation method is probably 
preferable since considerable labor is involved in getting 
the j . For higher order functions the orthogonal poly- 
nomial method will be necessary. The two-dimensional poly- 
nomial case will be discussed next and then spline functions 
will be introduced. 
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Two Dinvensional Polynomial Approximation 

Two dimensional approximation is a direct extension of 
the one dimensional case for the norm. The chief problem 
lies in the size of the problem in that the number of terms 
increases as the square of the order of the approximation. 

The two dimensional image distortions Ax and Ay must 
be estimated for every element in the picture from a limited 
number of unequally spaced points at which the distortion is 
known. The functional approximation problem is solved here 
by computing a least squares polynomial approximation to 
the given points using the normal equations. Given are a 
set of displacement values distributed over the two dimen- 
sional image space. Let f y^ x i' y i^' «* n be 

the true image displacements at the point x^,y^. Tbe displace- 
ment over the entire image is to be approximated by a polynomial 
function based on the n measured values: 


m m 

ax(x,y)= k I 0 a jk*jk (x ' y> 

{ 19 ) 

m m 

y x ' yl “ 3=0 iJo b jk' , ’jk (x ' y) 


For approximations of low order the nonorthogonal 
monomial basis function x-* can be used and the normal 
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equations solved for the desired coefficients. The 
imating functions in this case become: 

Ax(*,y)= a o 0 +a i 0 x + a ol y + a 20 x2 '" 

mm 4 
m y, r a vJ,, k 

jSo k-o a jk x - 


Ay(x,y)» b Q0 + b lft x + b ni y + b on X 


'10 


01 J 


' 20 ' 




mm . . 

■ j£o kJo b jk x y 

where: Ax,y are the approximated values 

{a,b} are parameters to be determined 
The elements of the equation matrix G become: 

1 f • . « n 

A-l,.. (m+ 1) 2 

j=0, . . . ,m 

k 3 Of. . . f m 

And the normal equations are again generated by 

T T 
G 1 GA«G A f x 

m m 

G^GB-G^fy 

where: A,B are column vectors of the desired 

coefficients a^, bjy. 

f x ,fy are column vectors of the image 
distortions in the x and y directions. 


appro?*-* 


( 20 ) 


( 21 ) 
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The solution for the coefficients of the least squares 
functions are 

A- (G T G) -1 G T f x (22 

B* (G T G) -1 G T f y 

The ordering of the powers of x and y in the polynomial and 
the elimination of certain terms are two items of variation 
in procedure. The coefficient vector is actually a matrix 
of terms 


a 00 a 01 a 0m 


a 10 a ll 


[ a m0 • • • a mm J 

The subscript notation is chosen so that the first subscript 
is the power of the x term and the second is the power of y. 
And the set of polynomials is also a matrix which in the 
case of <J> j (x) *x^ is 


“ 1 

x x 2 ... 

x m 

y 

2 2 

xy x y ... 

*v 

y 2 



_ m 


m m 


and the approximating function is expressed as 


Ax(x,y)« ♦ (x,y) §A 
Ay (x,y)= ♦(x,y)«B 


{ 23 ) 
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Where the (g) indicates the matrix dot product: 

m m 

A (x) B a ij ,b ij* The normal equations are 

solved in the same way as for the single dimensional caise 
except two coefficient sets are obtained, one for Ax and 
the other for Ay. 

Two Dimensional Orthogonal Polynomial Approximations 

The generation of two dimensional orthogonal polynomials 
is again a straight-forward extension of the development 
for the one dimensional case. For the case of unequally 
spaced points, a recursive relationship can be used to generate 
the polynomials. Polynomial can be obtained as a function 
of lower order polynomials: 

4 >j + 1 ( x 'y>= <x+y-dj +| )*j < x 'y) (24) 

where ^ (x,y)=Q and <J> (x,y)=l 
1 0 

we require that 
n 

i=!*j ( x i' y i> *k+i ( x i' 7 i ) =0 for D = 0 r i # • • • 


( 25 ) 
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Substituting the recursive expression into this requirement 
for j=k gives : 

n n 


i i, ♦j (x i'5'i )x i*k (x i'i'i ) ♦j< x i'Vi>y ♦k (x i' y i 


n 


n 


for j*0, 1, . . .k (26) 

A double summation over (x^,yj) is not used because each 
coordinate pair is unique due to the random distribution 
of x^ and y^ thus the single summation takes into account 
all possible points in the two dimensional space. 

For j>=0 , 1 , . . «k-2 in equation 26 the rightmost two terms 
are zero because the orthogonality condition holds for poly- 
nomials up to the one being generated. 

The first two summations with x^ and y^ in the 
summands are polynomials of degree not greater than k-1 and 
can be expressed as a linear combination of the (x,y) and 
thus these terms will be zero also. For the case j=k~l the 
third term is still zero but the fourth is non-zaro. The 
first two terms are non-zero because the order of 

x i^k-l * x i ' y i^ * and P roduct summation with $^(x,y) will 

be non-zero. Thus a condition for g is created: 

B k- i£ ♦k- J (x i'yi )x i*k (x i^i )+ iIi*k- 1 (x i'yi ) yi*k (x i' y i ) 


i«i [ *k-i (x i' y i } ]1 


( 21 ) 
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For j=k the fourth term is zero and the third term is non- 
zero. This gives an expression for a: 
n n 


l k+i' 


,E 

i=i 






( 20 ) 


n 

I 


W, 


(Xi^f) ] 


Thus a , can be found which generates a sequence of ortho- 
gonal polynomials. Given that 4^ = 0, = 1 

♦l (x,y) = (x+y-ot^) . 

It is required that 
n 

i | J ^ 0 (x,y)4> 1 (x,y)=0 

or (’‘ i + S' i -« 1 >=° 
n 

Thus ou =1 .£ (x.+y.) 
x n 1-1 1 

From this and are orthogonal and is generated from 
these two polynomials known to be orthogonal thus the 
assumptions used in deriving a and 6 will hold. The 
generation of the approximating function from the is 
identical to the case for one dimension. 
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m 


where: 


Ax(x,y)= jfio a j*j * x,y) 
m 

AY(x,y)«= ji 0 b j*j 

a j = i£> x iVV y i/i.VV y i 


( 29 ) 


) , similarly for b ^ 


Use of Tensor Products for Two Dimensional Approximations 
For a coordinate system in which the independent 
coordinates of the given date form a Cartesian product set 
X(3£>Y the following tensor product of functions is defined: 
The tensor product of two sets of functions 
{$^|i= 1,2,. ..,p} and {i^|i=* 1,2, ...q) is the set 

| l-i*P, l=j=q>. This product is stated as 

{*i> cb {*j>- 


The linear approximating function formed by such a 
product is: 

m m 

F(A,x,y)« ^ a i b j^i. <y) (3 ° 

Where: A represents the coefficients a,b 

m is the order of the function 
It can be shown that if the functions {4>j_ j i=l , • . . ,m} 
form an orthonormal set then the tensor product is an 
orthonormal set. Then the best L approximation (least 
squares) can be found by solving the normal equations using 
the tensor product form of function. 
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The tensor product approach eliminates the need for 
constructing the two dimensional orthogonal polynomials which 
becomes an extremely cumbersome task. The tensor approach is 
used in the spline function method to be discussed next. The 
spline functions contain as a subset the polynomial function 
approach discussed above; thus, the following sections cover 
both cases and represent the major thrust of this report. 

Spline Function Approximation 

The one and two-dimensional approximating functions dis- 
cussed up to this point are assumed to be valid over the entire 
region of interest with adequate accuracy. For cases in which 
the order of the distortion is low and the area covered is 
limited, the single function approach is adequate. When the 
image distortions become severe and the area to be represented 
increases the order of the functions required becomes imprac- 
tically high. Use of two-dimensional functions of fourth, fifth 
or higher order is undesirable for computational reasons and 
because the number of control points required becomes excessive. 

The bi-quadi.atic function is of the form: 

2 2 2 2 2 2 

Ax(x,y)= a +a x+a x+a xy+a x +a y +a x y+a xy +a x y ( 31 ) 

012 3 ** 3 * 7 • 

and requires nine coefficients. The bi-cubic function requires 
16 coefficients. Representation of higher order curves or sur- 
faces can be achieved through use of lower order polynomials 
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joined together in a piecewise continuous manner. Such approx- 
imations are called spline functions [6] and constitute a class 
of extremely useful and successful functions which were first 
considered from a mathematical viewpoint by Schoenberg in 
1946 [9] and research on these functions has increased steadily 
since. The great value of the piecewise polynomial approach is 
that complex disjointed functions arising from real, physical 
situations can be approximated rather conveniently. The random 
variations in multispectral scanner image geometry are generally 
unrelated from one place to another. Whereas for polynomials 
and most other mathematical functions their behavior in a small 
region determines their behavior everywhere. Piecewise contin- 
uous (spline) functions do not have this problem and for cubic 
or higher order polynomial splines the splines are smooth curves 
in the physical world. 

A one-dimensional spline function is defined by a set of 
points called knots 

*• < ^k < ^k+l“ b 

over the interval la,b] and a set of polynomials of degree n 
valid between the knots and having n-1 continuous derivatives 
at the knots. One representation as presented by Rice [6] of 

splines is of the form: 

k n 

S(A,E,x)= 2 a. (x-£.) +JI(x) 

i=l 1 1 + 


( 32 ) 


- 20 - 


Where : 



II (x) = Polynomial of degree n with coefficients 

a^,i=k+l, ... k+n+1 

A = Parameter Vector (a^,a ? , ... 

H = Set of knots ... 5^^) 

Man} other forms of representation exist for splines; however, 
the resulting functions are the same. Splines of greater them 
third order are generally not used, the advantage of the spline 
approach being that high order polynomials can be avoided. A 
first order spline would be a sequence of linear or ramp func- 
tions joined together at the knots forming the familiar piece- 
wise linear functions. A second order spline would be a set of 
quadratic polynomials connected at the knots and having contin- 
uous first derivatives at the knots. Similarly, a third order 
spline would have cubic polynomials joined at the knots and 
having continuous first and second derivatives. 

A more common representation for splines is: 

= Z C (x-5 ) 3 e -x-S 

j=0 ij i-1 i-1 i 


S(A,-,X) 


( 33 ) 
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which illustrates the piecewise polynomial nature more explic- 
itly. Also, when n is odd the splines can be uniquely repre- 
sented in terms of the values and derivatives at the knots. 
Thus, for linear and cubic splines, the approximation can be 
completely specified by the values: 



d j S(A,S,x) 



x-?i 


i=l, ... fkH, ^ = 0,1, 


. . . 


n-1 

T~ 


The 1 st and 3 rd (odd) order splines turn out to be extremely 
convenient choices since for even functions the number of 
derivatives needed to be specified at the left and right knot 
is not equal and the spline cannot be uniquely specified by 
knot point values and derivatives. Furthermore, since the linear 
splines have n-l=0 continuous derivatives at the knots the 
spline function has jumps in slope at the knots and do not form 
"smooth" curves. Thus, the cubic spline becomes the natural or 
preferred order for spline approximating functions. The cubic 
spline in one dimension is of the form: 


£(A, z ,x )= c ii +c i2 ^ x ~^i-l ^ +C i 3 ^ x- ^i-l ^ "^^4 (x— 1 ) 

for 


( 34 ) 


Note that four coefficients are required for each polynomial 
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which can thus be specified by four knot point values ; 


i-1 0 

S(A ^ E#5 i-l ) 

v i(0 = s(A.s,e i ) 


dS (A, - , x) 


dS (A, H , x) 

i-1,1 

dx 

x=5 i-i 

v i,r s 


x=£. 


It can also be shown [10} that the entire cubic spline function 
is uniquely specified by only the function values at the knots 
and derivatives only at the end point knots. 


Two Dimensional Spline Function s 

The discussion of splines thus far has been in the context 
of one-dimensi anal functions. Second and higher dimensional 
splines are more complex analogs of the one dimensional case. 

The two-dimensional case will be discussed in terms of cubic, 
or more precisely, bi-cubic polynomials since they have the 
same advantages that the one-dimensional cubic splines have. 

It can also be proven that specification of 16 corner values 
and derivatives uniquely specifies the spline function polynomials. 
The two dimensional splines are defined on a rectangular grid 
in a two-dimensional plane. Let the divisions between the 
rectangles or knots on the X axis be defined by the set {£^} 
and the knots on the Y axis by the set . Then the 2-dimen- 
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sional spline polynomials are of the form: 

3 3 

AX<x,y) * Z Z a^ (x-5 i „ 1 )* n {y-y, x ) n ( 35 ) 

m*0 n°0 J 

For £i_i“ x “£j_» yj_j^y^Mj i aa l»..« f K f j=l,...,L 

For each polynomial there are 16 defining coefficients. A 
spline grid is depicted in Figure 1 for K*3 and L*3. The 16 
corner conditions on each rectangle that specify the cubic 
polynomial follow Let the spline function be represented as 
S(x,y)j then the corner values for rectangle ' 


(5 i ,y j-l > ( 

' ^i-1' ; ' (5i»Mj) are: 



S^i-l'Hj-l* s(? i' M i. 

<!> 

Jj) S(q,y.) 


3S(x,y) 

3S (x,y ) 

3S(x,y) 

3S (x,y) 


3x 

x=q-i 5* 

x-5 £ 3x 

x=5 i-i 3x 

x =?i 


y “ M j-l 

y ” y j-l 


y=M : 

3S(x,y) 

35 (x,y) 

3S(x,y) 

3S (x,y) 


ay 

x=5 i-i ** 

5y ' 

x=5 i-i 3y 

x=5 i 


y =u m-l 

y=y j-l 

y =Mj 

!y=u j 

3S(x,y) 

3S (x f y) 

3S (x,y) 

3S (x,y) 


3x3y 

x-^i 1x5y 

x*^ 3x3y 

x*^ Sx3y 

x=C i 


y^j-i 

ycy j“l 

y=Wj 

yssy j 
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Fort unate ly the two dimensional analog of a previously stated 
result holds for the bi-cubic splines. That is that the 
polynomials are uniquely determined if only the value of the 
function is specified at each of the mesh points and the deriv- 
atives are specified only at the outside edges of the grid. 
Specifically, the piecewise bi-cubic function is uniquely 
specified by only the following values: 

i=0,...,L, ( (L+1MM+1) values) 

i=0,L, j=0,l,...,N (2M values) 

i=0,l,...L, j*0,M (2L values) 

i“0,L, j«0,M (4 values) 

This elegant result will not be proven here , but the proof can 
be found in the cited report by de Boor [10] . For the nine 
rectangle mesh depicted in Figure 1 specification of 16 mesh 
function values, 16 edge derivatives and four cross corner 
derivatives or 36 total quantities completely specify 


S ( t w j ) 

3S|x,y) 


x=^ 

y-M. 


3S(x,y) 

“^y 




3S(x,y) 


y=^ 
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the nine bi-cubic polynomials forming the spline function. The 
representation without the simplification would require the 16 
values, 16 X derivatives, 16 y derivatives, and 16 cross deriv- 
atives or 64 values. 

Computation of the bi-cubic polynomial coefficients from 
the specified mesh point values can be accomplished through 
simultaneous equation solution also described in [10], The 
method discussed assumes knowledge of the exact values and 
derivatives of the function at the grid nodes. In practice, 
these values are not known exactly and generally are a set of 
approximately known values and these are often unevenly spaced 
over the domain of the function. Thus, the real problem is 
computing an approximate two dimensional spline function based 
on this data. The least squares spline problem is thus the 
problem to be solved. 

Least Squares Two-Dimensional Spline Approximation 

The least squares solution for the two-dimensional spline 
approximation function is most advantageously computed using a 
set of orthogonal spline basis functions. Using this approach, 
once an orthogonal spline basis is obtained the approximation is 
easily computed using inner products. If is a complete 

orthonormal spline basis for the set of all spline functions 
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over the domain of interest, then 

i,j“0,...N (36) 

Where: 6 i; .=l if i=j and sero otherwise. 

The approximating spline function S is computed as 
n 

S* S <f,^> <l> 4 (37) 

i=0 1 1 

If a non-orthogonal basis { 4>£ J is known the orthonormal 
basis {ij^} can be computed using a procedure similar to the 
one described for the one dimensional case. The Gram-Schmidt 
orthonormalisation procedure is usually used for this purpose. 
This procedure is described by the two step iterative formula: 

c VIM 

i-l 

= 4>.- E <4>. ,iK> ip . i=»0 , . . .N (38) 

1 1 j=0 1 3 3 

h m VIM 

By this process, the non-orthonormal basis { > is converted to 
the orthonormal spline basis function set {i|k} which can be 
used directly for evaluating the approximating function co- 
efficients. It is pointed out in (7) that in forming the initial 
basis it is advantageous to construct each 4 k so as to have 

one more extremum than This tends to generate a "nearly" 

orthogonal initial basis which improves the accuracy of the 
resulting orthogonal basis. 
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The procedure for solving for the least squares bi-cubic 
approximating function consists of solving for the orthonormal 
spline basis functions for a given set of knots and then eval- 
uating the approximating function and the mean squared (L 2 ) 
error. The knots can then be moved or increased in number to 
attempt to reduce the error. An algorithm for carrying out this 
process if presented in [7] for approximation in one independent 
variable. Development of a two-dimensional cubic analog for 
this process was carried out as part of this project. The 
algorithm facilitates geometric correction and registration of 
aircraft scanner data and similar data having almost any degree 
of distortion. 

Two Dimensional Spline Function Approximation Algorithm 

The algorithm developed is a generalization to two dimensions 
of the algorithm FXDKNT described in [7] by de Boor and Rice. 

The tensor product approach is used in generation of the spline 
basis functions rather than attempt to compute bi-variate basis 
functions. The algorithm was originally written to accomodste 
100 data values to be approximated and up to 26 knots in addition 
to the left and right boundary knots, or a total of 28. In the 
two dimensional version the number of points was kept at 100 and 
the number of Y axis knots made the same as for X or 28. This 
greatly expands the size of the program but keeps the same capa- 
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bility in the new program on each axis should it be needed. In 
practice the number of knots on each axis for scanner imagery 
probably will not exceed five or six; however, in other applica- 
tions the full power of the algorithm may prove useful. 

The tensor product form of basis function generation results 
in two sets of orthonormal spline basis functions (i^(x)} and 
{}i j (y) } from which the coefficients are obtained for the ortho- 
gonal projection of the function to be • approximated onto these 
basis functions. The form of the approximating function is thus: 
IX IY 

AX(x,y)« Z £ a. (x) y . (y) (39) 

i-lj-1 13 1 3 

Where ere the spline basis functions 

the coefficients of the approximation 

function in the basis i p. (x) y. (y) 

1 1 

IX, I Y are the number of X and Y dimension basis 
functions in the solution. 

The are computed as the inner product of the function 
to be approximated and the basis functions: 

a i j= <f x (x,y) V' i (x)Uj(y)> (4°) 

Where 

< > denotes the two dimensional inner product 
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The algorithm proceeds by first computing a single cubic 
polynomial approximating function over the entire set of points 
to be approximated. Then additional knots are introduced one 
by one and additional basis functions are computed. The coor- 
dinate of the new basis function is computed and the contribu- 
tion of the new term is subtracted from the remaining error in 
the approximation . 

The polynomial spline functions are computed from the basis 
functions and their coordinates by evaluating the value , deriva- 
tives in the x,y and cross directions at the corners of each 
of the rectangular segments of the domain being covered. These 
sixteen values are then used to define a cubic polynomial for 
each rectangular region. The sixteen values are then transformed 
into the 16 polynomial coefficients. The economizing procedure 
discussed previously is not used in the current algorithm. The 
resulting approximation is represented by the sixteen coefficients 
for each spline region for as many regions as were specified by 
the knot set. Thus, a function having four knots in the x 
direction and six knots in the y direction, including boundary 
knots, would have (4-1) x (6-1) xl6 = 240 quantities specifying 
the approximating function plus the ten knot values . 

The algorithm can be re-executed to add or delete knots to 
adjust the overall RMS error in the approximation to a desired 
level. An artifice was used in the computation of the two 
dimensional inner product to handle the case of randomly located 
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data points. The data points are constrained to lie on a 
quasi-re ctangular grid and the means of the resulting groups 
are used as the x and y abscissa values in the inner product. 
A nearest neighbor rule is used to assign function values to 
the points at the nodes of the artificial uniform grid. This 
enables a simple trapezoidal integration inner product to be 
computed but causes error in the approximation. An iterative 
technique is then employed to correct for this error. 
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Example of Spline Function Approximation 

The multispectral aircraft scanner system flown by the 
ERIM* organization produces imagery in long strips of nominally 
two miles in width at 5,000 feet altitude. This data is often, 
subject to severe distortions due to pitch and yaw variations 
in the aircraft attitude and lateral motions due to cross winds 
since the scanner is fixed to the frame of the aircraft. The 
scanner is roll stabilized so that only pitch and yaw angular 
distortions are experienced. Thus, this type of imagery can 
be affected by five platform variables: pitch, yaw, and trans- 

lational variations in three dimensions. An example of air- 
craft motion distortion in the MSS imagery is shown in Figure 2. 
More sophisticated scanners are stabilized on the pitch, roll, 
and yaw axes; however, this requires costly gimballing mounts 
and costly support control systems. In most scanner imagery 
cases, some degree of random distortion will be present and the 
spline function techniques are expected to be useful in a wide 
range of cases. The extreme flexibility of the spline approx- 
imating functions allows the case of using only one function 
for the entire image for simple distortion up to the case of 
many spline function regions covering the image to be handled 
with the same algorithm semi-automatically . 


* Environmental Research Institute of Michigan 
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Fiqure 2. Aircraft Multi- 
spectral Scanner iraaery 
(.5R-.65U band) showina 
severe geonetric distortion 
due to crosswinds. LAPS 
Pun No. 71054100. Area is 
imnediately west of Craw- 
fordsville, Indiana. 

Date: August 17, 1971. 

Altitude: 5,000 ft. 




Fiqure 3. Aerial photo 
of area in Fioure 2 
showinn the desired 
imaqe cteoiretry. 
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The multispectral scanner imagery shown in figure 2 in 
a typical example of EUIM low altitude flight data. The aerial 
photograph segment shov:n in Figure 3 covers the area imaged by 
the scanner and represents the desired geometric shape of image 
in Figure 2. The function required to transform the MPB imagery 
into the geometric form of the photo is specified by defining 
checkpoints or matching points in the image and map. These 
points can be obtained by a variety of manual or automatic 
methods and for this example they were obtained by measuring 
the coordinates in inches on the MSS image and photo using a 
coordinate digitizing table. The scale of the image is approx- 
imately 1:56300 and the scale of the photo is approximately 
the same. The coordinates of the checkpoints digitized from 
the imagery and photo are listed in Table 1. 

The values from Table I were input to the two dimensional 
cubic spline algorithm first for the case of only one block. 

This results in a cubic polynomial fit over the entire region 
which in this example was for .75 <. x < 2.594 and .593 < y < 
10.25. The results of the approximation are listed in Table 2 
which includes the values to be approximated (the x and y posi- 
tion of each conjugate point in the aircraft scanner image) , 
the approximations, the error and three error statistics. The 
root of the mean squared error for the approximation is .063 
for the x coordinate and .11 for y. The maximum error was .132 
for x and .37 for y. 
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Table 1. Coordinates of Matching Points for Aircraft 
Scanner Data Correction Example. Purdue 
Flight Line 212. Scanner data obtained 
august 17, 1971. Aerial photograph made 
from Color IR photograph taken at 60,000 
feet by NASA RB-57 in 1971. 


Point 

No. 

Photograph 

Scanner 

Image 

X 

Y 

X 

Y 

1 

1.53 

.594 

1.06 

.75 

2 

.781 

.813 

.063 

. 875 

3 

2.125 

.593 

1.875 

.813 

4 

2.438 

1.188 

2.125 

1.375 

5 

1.563 

1.188 

1.0 

1.031 

6 

.750 

1.188 

.031 

1.188 

7 

.969 

2.188 

.188 

2.250 

8 

1.563 

2.250 

. 875 

2.250 

9 

2.250 

2.656 

1.875 

2.625 

10 

1.031 

4.093 

.438 

3.938 

11 

1.562 

4.813 

1.0 

4.063 

12 

2.125 

4.25 

1.781 

4.094 

13 

,969 

5.375 

.25 

5.188 

14 

1.531 

5.375 

1.031 

5.188 

15 

2.281 

5.344 

1.969 

5.250 

16 

.968 

6.50 

.188 

6.313 

17 

1.565 

6.50 

.938 

6.313 

18 

2.438 

6.50 

2.0 

6.375 

19 

.938 

7.656 

.125 

7.469 

20 

1.531 

7.688 

.813 

7.50 

21 

2.50 

7.656 

2.063 

7.50 

22 

1.094 

10.25 

.156 

10.031 

23 

1.625 

9.969 

.75 

9.75 

24 

2.375 

9.938 

2.125 

9.75 

25 

2.031 

.562 

1.75 

.781 

26 

2.0 

1.125 

1.70 

1.312 

27 

1.875 

2.5 

1.375 

2.47 

28 

2.125 

4.75 

1.81 

4.65 

29 

1.781 

5.375 

1.437 

5.218 

30 

2.125 

6.5 

1.70 

6.375 

31 

2.125 

7.656 

1.65 

7.5 

32 

1.875 

9.94 

1.187 

9.75 


• 

ft 
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Table 2a Error data for x dimension approximation 

usinq one cubic snline olock^ 

KOIJT MILAN SGUARE ERROR = 0.6257/9E-G1 

AVERAGE ERROR = 0.52AU7HE-01 

MAXIMUM ERKCJR = G.132832E X Ai 1.780999 5.37A999 

APPROXIMATION AND bCALEI) ERROR CURVE 

DATA POINT APPROXIMATION UFVIATIUN X lCt + 2 

1 

0.7A99997U 

1 . 18799973 

0.03031971 

-O.C6U076 

2 

0 . 78099966 

C. 8129996/ 

-0.U66A 7876 

-12.9A7B15 

J 

0.93799961 

7.65599918 

0. IMWH 

1.27A5A9 

A 

0. 9679995b 

6. A9999905 

0.2262 7 86 J 

3.827875 

5 

0. 96R99956 

2 . 1879997 j 

0.26086875 

7. 2868A0 

6 

0.96899956 

5. 37A9990 1» 

0.2/540171 

2.5A0183 

7 

1.03099918 

A.092999A6 

0. 33013874 

-10. 786 1C A 

0 

1.09399986 

10. 2A999906 

0 . 0 9 0 2 2 1 A 1 

-6.57 78 A9 

9 

1.52999973 

U . 5939y96d 

0 . 9A78 3 A 7 J 

-11.216A7A 

10 

1 .53099918 

5 . 37A99905 

0.91989422 

-11. 110497 

11 

1.53099918 

7.68799877 

0.86078262 

A. 778296 

12 

1.66199932 

A. 61299877 

0.9659261 7 

-3. A 07383 

* 13 

1.56299973 

1 .18799973 

0.97424936 

-2 . 5 7 506 A 

1 A 

1.56299973 

2 . 2 A999905 

C. 95527613 

8.027655 

15 

1 .56A99950 

6. A9999905 
9.96899891 

0.95872116 

2 . G 7 2 1 5 5 

16 

1.62 A 9 9905 

0 . 7 A 2 0 5 3 9 9 

-0 .794572 

17 

1.70099918 

5. 37A99905 

1.30 A 16775 

-13. 2B 3157 

1H 

1.8 7 A99905 

2 . A9999905 

l . 4 58 5 L 8° 3 

8.351898 

19 

1.8 7A 99905 

9.93999958 

1 . 1A7&89*2 

-3.930950 

20 

2.00000000 

1. 12A99905 

1.64196968 

-5.802917 

21 

2.03099918 

0.561 9 y97A 

1.68718338 

-6.281567 

22 

2. 12A99905 

A . 7A999905 

1.77307/96 

-3.692150 

23 

2. 12A99905 

6. A9999905 

1.73779392 

3.779507 

2A 

2. 12A99905 

7.6559991b 

I . 6965A8Ao 

A.65A88A 

25 

2 . 12 A99905 

0.592999 70' 

1 . 81 5bti 1 A o 

-5.931759 

26 

2 . 12AV9905 

A . 2A9 99905 

1.77913666 

-C .186253 

27 

2.2A69990D 

2.65599918 

1.90832806 

3.3329C1 

2d 

2.28099918 

5. 34399891 

1.90369275 

-6.530666 

29 

2.37A99905 

9. 93799877 

1.912562A j 

-2.54163 7 

30 

2. A 3 799973 

1.18799973 

2. 1 1 A 6 5 5 A 9 

-1 .034355 

31 

2.A37-99973 

6. A 9 999905 

2 . 0 1200867 

1.200867 

32 

ft- 

2. A9999905 

7.65599916 

2.09133911 

2.833939 
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Table 2b Error data for y dimension approximation 
using one cubic spline block. 


ROUT MEAN SQUARE ERROR = 
AVERAGE ERROR 
MAXIMUM ERROR = 


C.llObilE 00 
0.76426 1 C-0 1 

0.360560E 00 A I l .6ol444 


APPROX I MA n ON AND oCALED ERROR CUR VE 

DATA POINT APPROXIMATION DEVIATION X 

QuCiOTA 1 . I fl7au<17 4 I . Xjy 4 UT1>\ I A.hlx'ti 


1 

0. 74999970 

1 . 187999 73 

1.32238770 

13.4 3b7 I 7 

2 

0.78099966 

0.81299967 

0.65502316 

-21 .99 76 i b 

3 

0.93799961 

7.6559991a 

7.442 79766 

-2.620125 

4 

0.96799958 

6.49999909 

6 . 322896° 3 

0.9896.9 

3 

0.96899956 

2.18799973 

2.182624*2 

-6.7 3 7* 2 3 

6 

0.96899956 

5.37499905 

5.28064754 

6.264/02 

7 

1 .03099918 

4 .09299940 

3.88046741 

-5. 753136 

8 

1.09399986 

10.24999905 

10 .064022°6 

3.3*2288 

9 

1.52999973 

0.5939996O 

0.80707961 

5. 707911 

10 

1.53099918 

5.37499905 

5.002182°l 

-10.681772 

11 

1.53099918 

7.6879987/ 

7.49032593 

-0.96 7312 

12 

1.56199932 

4.81299877 

4.43155956 

36.8559/2 

13 

1 . 56299973 

1. 1879997J 

1.19262123 

16.162201 

14 

1.56299973 

2.24999905 

2.004449*4 

-24.654916 

15 

1.56499958 

6.49999905 

6« 2231597 i 

-8.903994 

16 

1.62499905 

9.96899B9t 

9.76933289 

L .933304 

17 

1 . 78099918 

5.3749990b 

5. 16595454 

-5.2 0 44 d 7 

18 

1.87499905 

2.49999906 

2.481965°/ 

1. 19 6570 

19 

} .87499905 

9.9399995b 

9./438964y 

— v. 6 1 620 6 

20 

2.00000000 

1.12499905 

1.2920284 i 

-1 .997089 

21 

2.03099918 

0.5619997* 

0. 79676783 

1.5/081 7 

22 

2.12499905 

4.74999905 

4.67579651 

2.579689 

23 

2 .12499905 

6.49999905 

6.378982*4 

0. 398300 

24 

2. 12499905 

7.6559991b 

7 • 5 1 420 5° 3 

L .4 2068 9 

25 

2.12499905 

0.5 92 9997^ 

0. 86055422 

4.765496 

26 

2.12499905 

4.24999905 

4. 1960 1 44 3 

10. 201404 

27 

2.24999905 

2.6559991b 

2.61439323 

- L ... 0650 L 

28 

2.28099918 

5.34399891 

5.16276550 

-8 . 72 3364 

29 

2.37499905 

9.93799877 

9. 734436°* 

-l.osoin 

30 

2.43799973 

l . 18799973 

1.356300*5 

- 1.869804 

31 

2 .43799973 

6.49999905 

6.25756836 

-1 1 . 74306 9 

32 

2.49999905 

7.6559991b 

7.57215851 

7.2169/7 


*4 • ^ 


lul *-2 


(Jcxt the y dimension was divided at the approximate mid- 
point by an additional knot at y * 5.0 and the spline approx- 
imation was computed for the two regions. Two cubic spline 
functions were thus computed which join with continuity in 
value and first and second derivative at the line y =* 5.0. 

The fit was improved to an r.m.s. error of .04 for x but the y 
error remained about the same at .105, The maximum error was .008 
for x and .287 for y. The two section spline improved results 
considerably and produced a smooth curve with no discontinuity at 
the knot line. The results are tabulated in Table 3 listing the 
same information as Table 2. This is a simple illustrative ex- 
ample and no attempt will be made here to optimize the fit to 
the aircraft data by varying the position of the knot or adding 
more knots. An algorithm which optimizes the positions of the 
knots is being developed as a continuation of this work. A 
great deal of flexibility is available in the spline approach and 
the error could be further reduced by appropriate manipulation 
of the knots. 

Summary and Conclusions 

This report presents a discussion of least squares approx- 
imation techniques with two dimensional spline function approx- 
imation being the main topic. A one dimensional algorithm 
due to de Boor and Rice was described and its extension to two 
dimensions is the subject of the work reported here. The algo- 
rithm is operational; however, certain problems with the two 
dimensional inner product remain to be solved. A technique 
was used in the current program in which a nearest neighbor 
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Table 3a Error data for x dimension approximation 
using two cubic spline blocks with new y 
knot at 5.0. 


ROOT MU AN 5WUARE ERROR = 

average error = 

MAXIMUM ERROR 


0.404584E-01 
0. 346769E-01 

C .888879E-01 AT l .t 1099 / 4 .<,92 ' r 


1 

APPROXIMATION AND bCALED 
DATA POINT 

0.74999970 1.18799973 

ERROR CURVE 

approximation 

0.03919234 

DEVIATION X 
8. 192368 

2 

0.7809996b 

0.8129996/ 

0.01602793 

-46.971954 

3 

0.93799961 

7 ,65599918 

0.09274996 

-32 . <i 5U 6 j 9 

4 

0 . 96799958 

6.49999905 

0. 24594599 

67.946136 

5 

0.9689995b 

2.18799973 

0.210124^3 

22.124222 

6 

0.96899956 

5.37499905 

0,33351320 

83.613356 

7 

1 .03099918 

4.09299946 

0.34911191 

-68.887863 

8 

I.C9399986 

iC. 24999903 

0.11023430 

-46. 765666 

9 

1.52999973 

0.59399968 

1.09113693 

31.137406 

10 

1.33099918 

5. 37499903 

0.97816896 

-62.830215 

11 

1 .53099918 

7.6879987/ 

0.80351156 

-9,466106 

12 

1.56199932 

4.81299877 

1.02039337 

20.393372 

13 

1.56299973 

1.18799973 

0.985126 7 1> 

- 14.97 3/4 3 

14 

1.56299973 

2.24999903 

0.87803519 

3 • u 3 562 4 

15 

l.! 36499958 

6.49999903 

0.972326^4 

34.32643 1 

16 

1.62499905 

9.96899891 

0.77036572 

20 .3660 i 3 

17 

1. 78099918 

5.37499906 

1.3768343/ 

— 60 . 16 444 4 

18 

1.87499905 

2.49999905 

1.37457943 

-0.4 1961 / 

19 

1.87499905 

9.93999956 

1. 16 14370 i 

-25.562266 

20 

2 . CO00Q000 

l. 1249990; 

1.657094 r O 

-42.904846 

21 

2.03099918 

0.56199974 

1.82848072 

/ 6 . 4 6 L 6 7 4 

22 

2.12499903 

4. 7499990; 

1.62524395 

15.244464 

23 

2.12499903 

6.49999905 

1.747889^2 

47.690656 

24 

2.12499905 

7.65599918 

1.639124P / 

-10.67474 8 

25 

2.12499905 

0. 592999/0 

1.9199343/ 

44.935226 

26 

2.12499903 

4.24999903 

1.80994034 

28.94114 / 

27 

2.24999905 

2.65599918 

i.b 7 827873 

3.279636 

28 

2.28099918 

5.34399891 

1.927209A5 

-41 . 769C4 7 

29 

2. 37499905 

9.9379987/ 

1.91138649 

-26.6122/4 

30 

2.43799973 
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20. 60 7266 
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Table 3b Error data for y dimension approximation 
using two cubic spline blocks witli new y 
knot at 5.0. 
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rule is used to define the values of randomly spaced data 
points at the nodes of a uniform grid. This makes computing 
the integral for the inner product simple but results in error 
in the approximation. A rule must be employed v?hen using this 
program in selecting data points over the two dimensional sur- 
face so that a quasi-uniform grid is maintained. The points are 
then grouped by the program and the mean of the group on the x 
or y axes is taken as the respective abscissa. Further work 
needs to be done on this and other problems relating to randomly 
spaced points in two dimensional approximation problems. Sub- 
sequent reports will document the algorithm in detail and address 
certain of these problems. 

It must be pointed out that the spline function approach to 
approximation is only one of a large number of methods each with 
their own advantages and disadvantages. For the problem of 
multidimensional approximation of functions the Weighting Function 
Technique of Jancaitis and Junkens [11] bears particular note 
and future work will evaluate this and other methods relative 
to the spline function approach. 

Finally, it should be noted that the multidimensional approx- 
imation techniques have application in many earth resources data 
processing areas in addition to image geometric distortion 
representation. Any case in which randomly located measurements 
are made of physical, electromagnetic, socio-economic processes 
is a candidate for this type of approximation technique. Speci- 
fically, the conversion of digitized topographic contours, 


airborne radiometer, and magnetometer and other geophysical 
data to a uniform grid format for computer overlay and image 
processing and analysis is the next application goal of the 
present work. Progress in these areas will be reported in 
subsequent documents. 
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